"""
Scalar field implementation.
"""

using LinearAlgebra
using StaticArrays
using Statistics
using ..Grids

"""
    ScalarField{T<:Real} <: AbstractDataField

Scalar field on a grid containing real-valued data.

# Fields
- `_grid::AbstractGrid`: The underlying spatial grid
- `_data::AbstractArray{T}`: The field data array

# Examples
```julia
# Create from function
grid = UnitGrid([64, 64])
field = ScalarField(x -> exp(-(x[1]^2 + x[2]^2)), grid)

# Create from data array
data = rand(64, 64)
field = ScalarField(data, grid)

# Create with zeros
field = ScalarField(Float64, grid)
```
"""
struct ScalarField{T<:Real} <: AbstractDataField
    _grid::AbstractGrid
    _data::AbstractArray{T}
    
    function ScalarField{T}(grid::AbstractGrid, data::AbstractArray{T}) where {T<:Real}
        size(data) == size(grid) || throw(DimensionMismatch("Data size must match grid size"))
        new{T}(grid, data)
    end
end

# Constructors
"""
    ScalarField(data::AbstractArray{T}, grid::AbstractGrid) where {T<:Real}

Create scalar field from data array and grid.
"""
function ScalarField(data::AbstractArray{T}, grid::AbstractGrid) where {T<:Real}
    ScalarField{T}(grid, copy(data))
end

"""
    ScalarField(func::Function, grid::AbstractGrid, ::Type{T}=Float64) where {T<:Real}

Create scalar field by evaluating function at grid points.
"""
function ScalarField(func::Function, grid::AbstractGrid, ::Type{T}=Float64) where {T<:Real}
    data = Array{T}(undef, size(grid))
    
    for idx in CartesianIndices(size(grid))
        coord = coordinate(grid, Tuple(idx))
        data[idx] = T(func(coord))
    end
    
    ScalarField{T}(grid, data)
end

"""
    ScalarField(::Type{T}, grid::AbstractGrid) where {T<:Real}

Create zero-initialized scalar field of given type.
"""
function ScalarField(::Type{T}, grid::AbstractGrid) where {T<:Real}
    data = zeros(T, size(grid))
    ScalarField{T}(grid, data)
end

"""
    ScalarField(grid::AbstractGrid)

Create zero-initialized scalar field with Float64 type.
"""
ScalarField(grid::AbstractGrid) = ScalarField(Float64, grid)

# Interface implementation
grid(field::ScalarField) = field._grid
data(field::ScalarField) = field._data

function Base.similar(field::ScalarField{T}) where T
    ScalarField{T}(field._grid, similar(field._data))
end

function Base.copy(field::ScalarField{T}) where T
    ScalarField{T}(field._grid, copy(field._data))
end

# Scalar-specific operations
"""
    gradient(field::ScalarField) -> VectorField

Compute the gradient of a scalar field.
"""
function gradient(field::ScalarField{T}) where T
    g = grid(field)
    dims = ndim(g)
    
    # Create vector field to store gradient
    grad_data = Array{T}(undef, size(g)..., dims)
    
    for d in 1:dims
        grad_data[ntuple(i -> i == dims + 1 ? d : Colon(), dims + 1)...] = 
            grid_gradient(g, data(field), d)
    end
    
    VectorField{T}(g, grad_data)
end

"""
    laplacian(field::ScalarField) -> ScalarField

Compute the Laplacian of a scalar field.
"""
function laplacian(field::ScalarField{T}) where T
    g = grid(field)
    lap_data = grid_laplace(g, data(field))
    ScalarField{T}(g, lap_data)
end

"""
    integrate(field::ScalarField) -> Float64

Integrate the scalar field over the entire domain.
"""
function integrate(field::ScalarField)
    g = grid(field)
    
    # For specialized grids with volume weighting
    if isa(g, PolarSymGrid) || isa(g, SphericalSymGrid)
        return integrate(g, data(field))
    elseif isa(g, CylindricalSymGrid)
        return integrate(g, data(field))
    else
        # Standard Cartesian integration
        vol_element = prod(spacing(g))
        return sum(data(field)) * vol_element
    end
end

"""
    get_extrema(field::ScalarField) -> Tuple{Float64, Float64}

Get minimum and maximum values of the field.
"""
function get_extrema(field::ScalarField)
    return extrema(data(field))
end

"""
    get_statistics(field::ScalarField) -> NamedTuple

Get statistical information about the field.
"""
function get_statistics(field::ScalarField)
    d = data(field)
    return (
        min = minimum(d),
        max = maximum(d),
        mean = mean(d),
        std = std(d),
        var = var(d),
        integral = integrate(field),
        norm = norm(field)
    )
end

"""
    interpolate(field::ScalarField, points) -> Vector

Interpolate field values at given points.
"""
function interpolate(field::ScalarField{T}, points) where T
    g = grid(field)
    values = Vector{T}(undef, length(points))
    
    for (i, point) in enumerate(points)
        # Find surrounding grid points and interpolate
        # This is a simplified implementation - full version would use proper interpolation
        idx = index(g, point)
        
        # Clamp to grid bounds
        clamped_idx = ntuple(ndim(g)) do d
            clamp(idx[d], 1, size(g)[d])
        end
        
        values[i] = data(field)[clamped_idx...]
    end
    
    return values
end

"""
    sample_at(field::ScalarField, point) -> Float64

Sample field value at a single point with interpolation.
"""
function sample_at(field::ScalarField, point)
    return interpolate(field, [point])[1]
end

"""
    get_slice(field::ScalarField, dim::Int, index::Int) -> Vector

Extract a slice of the field along given dimension at specified index.
"""
function get_slice(field::ScalarField, dim::Int, index::Int)
    dims = ndim(grid(field))
    indices = ntuple(dims) do d
        d == dim ? index : Colon()
    end
    return data(field)[indices...]
end

"""
    get_profile(field::ScalarField, start_point, end_point, num_points=100) -> Tuple{Vector, Vector}

Extract field values along a line from start_point to end_point.
"""
function get_profile(field::ScalarField, start_point, end_point, num_points=100)
    points = [start_point .+ t .* (end_point .- start_point) 
              for t in range(0, 1, length=num_points)]
    
    distances = [norm(p .- start_point) for p in points]
    values = interpolate(field, points)
    
    return distances, values
end

"""
    find_peaks(field::ScalarField; threshold=nothing, min_distance=1) -> Vector{Tuple}

Find local maxima in the scalar field.
"""
function find_peaks(field::ScalarField; threshold=nothing, min_distance=1)
    d = data(field)
    g = grid(field)
    peaks = Tuple{Float64, CartesianIndex}[]
    
    # Simple peak finding - can be improved with proper algorithms
    for idx in CartesianIndices(size(d))
        value = d[idx]
        
        # Check threshold
        if threshold !== nothing && value < threshold
            continue
        end
        
        # Check if it's a local maximum
        is_peak = true
        for offset in CartesianIndices(ntuple(_ -> -1:1, ndims(d)))
            neighbor_idx = idx + offset
            
            # Skip center point
            if offset == CartesianIndex(ntuple(_ -> 0, ndims(d)))
                continue
            end
            
            # Check bounds
            if checkbounds(Bool, d, neighbor_idx) && d[neighbor_idx] >= value
                is_peak = false
                break
            end
        end
        
        if is_peak
            coord = coordinate(g, Tuple(idx))
            push!(peaks, (value, idx))
        end
    end
    
    return peaks
end

# Convenience methods for mathematical functions
"""
    exp(field::ScalarField) -> ScalarField

Element-wise exponential of scalar field.
"""
Base.exp(field::ScalarField{T}) where T = ScalarField{T}(grid(field), exp.(data(field)))

"""
    log(field::ScalarField) -> ScalarField

Element-wise natural logarithm of scalar field.
"""
Base.log(field::ScalarField{T}) where T = ScalarField{T}(grid(field), log.(data(field)))

"""
    sin(field::ScalarField) -> ScalarField

Element-wise sine of scalar field.
"""
Base.sin(field::ScalarField{T}) where T = ScalarField{T}(grid(field), sin.(data(field)))

"""
    cos(field::ScalarField) -> ScalarField

Element-wise cosine of scalar field.
"""
Base.cos(field::ScalarField{T}) where T = ScalarField{T}(grid(field), cos.(data(field)))

"""
    abs(field::ScalarField) -> ScalarField

Element-wise absolute value of scalar field.
"""
Base.abs(field::ScalarField{T}) where T = ScalarField{T}(grid(field), abs.(data(field)))

"""
    sqrt(field::ScalarField) -> ScalarField

Element-wise square root of scalar field.
"""
Base.sqrt(field::ScalarField{T}) where T = ScalarField{T}(grid(field), sqrt.(data(field)))

# Comparison operations that return ScalarField
"""
    >(field::ScalarField, value::Real) -> ScalarField{Bool}

Element-wise greater than comparison.
"""
function Base.:>(field::ScalarField, value::Real)
    ScalarField{Bool}(grid(field), data(field) .> value)
end

"""
    <(field::ScalarField, value::Real) -> ScalarField{Bool}

Element-wise less than comparison.
"""
function Base.:<(field::ScalarField, value::Real)
    ScalarField{Bool}(grid(field), data(field) .< value)
end